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Abstract 

One of the long standing challenges in molecular simulation is the description of interfaces. On the molecular 
length scale, finite size effects significantly influence the properties of the interface such as its interfacial ten- 
sion, which can be reliably investigated by molecular dynamics simulation of planar vapor-liquid interfaces. 
For the Lennard- Jones fluid, finite size effects are examined here by varying the thickness of the liquid slab. 
It is found that the surface tension and density in the center of the liquid region decreases significantly for 
thin slabs. The influence of the slab thickness on both the liquid density and the surface tension is found to 
scale with in terms of the slab thickness S, and a linear correlation between both effects is obtained. 

The results corroborate the analysis of Malijevsky and Jackson, J. Phys.: Cond. Mat. 24: 464121 (2012), 
who recently detected an analogous effect for the surface tension of liquid nanodroplets. 

Keywords: Surface tension, long range correction, finite size effects 



1. Introduction 

Molecular simulation is a well-established approach 
for the analysis of fiuid interfaces and their molecu- 
lar structure. Much work has been dedicated to the 
interfacial tension.—"— For a fluid interface, there are 
(at least) three different aspects in which its size can 
be varied, each of which may affect the interfacial 
tension: 

• curvature effects, depending on the local char- 
acteristic radii of curvature 

• capillary wave effects, depending on the range 
of wavelengths permitted by the morphology 
and size of the interface 

• confinement effects, which arise due to spatial 
restrictions imposed on a fiuid phase by one or 
several interfaces or walls 

According to the Tolman— approach, the interfacial 
tension of a nanodroplet deviates from that of a pla- 
nar interface due to its extremely curved shape.—"— 
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However, it should be noted that all three phenom- 
ena are present when the size of a droplet is varied: 
Smaller droplets have a higher curvature, a smaller 
range of capillary wavelengths, and a more signifi- 
cant deviation from bulk-like behaviour due to con- 
finement. In addition to curvature, the other effects 
might therefore also have a significant influence on 
the formation of droplets in a supersaturated va- 
por, where the size of the critical nucleus and the 
nucleation rate are of major interest.—"— A similar 
case is cavitation, where bubbles emerge in a liquid 
phase.i&-iS 

Spherically curved interfaces of droplets were 
simulated for the first time by molecular simu- 
lation in the early 1970s. Nonetheless, while 
curvature effects are relatively well-studied, there 
are only astonishingly few investigations related 
to the other size effects, which are also present 
in the case of planar interfaces. Several previous 
works address the influence of small simulation vol- 
umes,—"— which is usually discussed in terms of 
capillary wave effects.—"— The present study con- 
siders the influence of the liquid slab thickness, i.e. 
of confinement by two parallel planar vapor-liquid 
interfaces which are close to each other. This effect 
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was previously investigated by Weng et al.«22 who, 
however, did not find a systematic correlation. 

The computation of interfacial properties is al- 
ways done in a single simulation volume containing 
both phases, the liquid and the vapor phase, sep- 
arated by the interface. For the calculation of the 
bulk properties there are many other methods, like 
Grand Equilibrium,— NpT plus test particle,— or 
the Gibbs ensemble method.— The surface tension 
can be computed for example via the virial route or 
the surface free energy.— i^^"— The virial route is di- 
rectly related to the common approach for calculat- 
ing the pressure in a molecular simulation. 22,23^38,39 
It is known that the pressure in dependence of the 
density exhibits a van der Waals loop in the two 
phase region.^2iiL 

At interfaces, the long range contribution to 
the interaction potential plays an important role 
for all thermodynamic properties.—"— Nonetheless, 
there are also simulations applying a truncated and 
shifted potential, which neglects the whole long 
range contribution.—!^"— When dealing with a 
homogeneous system, long range corrections are 
only needed for energy and pressure, ^^'^^ while in 
an inhomogeneous configuration, also the dynam- 
ics of the systems needs to be appropriately cor- 
rected.—"— 
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Figure 2: Surface tension 7 over the temperature for large 
liquid slabs. Comparison of results of different authors: (red 
squares) This work, (blue triangles) Janecek, (red dia- 
monds) Holcomb et al.,— (brown circles) Potoff and Pana- 
giotopoulos, (black circles) in't Veld et al.,^ (blue stars), 
as well as Lopez Lcmus and Alejandre,— ; (solid line) regres- 
sion from Eq. Bt . 



In the present work, the influence of the liquid 
slab thickness on thermodynamic properties is dis- 
cussed. A suitable long range correction is used to 
obtain accurate and validated results. 




Figure 1: Dependence of the computed surface tension 7 on 
the employed cutoff radius rc, for the truncated and shifted 
LJ potential without long-range effects (LJTS) and for the 
full LJ potential with a long-range correction (LRC) accord- 
ing to Janecek, The simulations were carried out in the 
canonical ensemble with N = 2 048 particles at T = 0.7 us- 
ing a simulation box with an elongation of £y = 66 in the 
direction perpendicular to the vapor-liquid interfaces. 



2. Simulation method 

In this study, the Lennard- Jones potential (LJ) 



/ ^ \ 12 








\nj) 





(1) 



is employed, where e and a are the energy and size 
parameters and is the distance between the two 
particles i and j. The standard LJ parameters e — 1 
and a = 1 are used, yielding a reduced LJ potential. 

As usual, the potential was truncated in order 
to reduce the computing time. A cutoff radius of 
Tc = 3 was used for the present simulations. To di- 
minish the error made by this assumption, a bin 
based tail correction was applied to the simula- 
tion.-^P'-- Thereby, the potential energy, the forces 
acting on the molecules, and the virial are each split 
into an explicitly computed part and a long range 
correction. The calculation of the correction terms 
was conducted every 10 time steps. It is known 
that this method provides cutoff independent re- 
sults for the LJ fluid. For a discussion of the em- 
ployed method in full detail, the reader is referred to 
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Figure 3: Density p (top) and differential pressure pjv — PT 
(bottom) over the y coordinate (i.e. the direction perpendic- 
ular to the interface). The temperature is T = 0.7. The 
blue dashed line corresponds to the minimum stable config- 
uration which is S = 4.3 for this temperature, while the red 
dash dotted one corresponds to 5 = 7 and the black solid 
one to S = 12. The dotted line in the upper picture rep- 
resents the bulk liquid density and the difference between 
the vertical dashed lines in the upper picture represent the 
equimolar slab thickness. 



Janecek's work,.^ wherein this approach was first 
presented. 

The interfacial tension 7 is given by the differ- 
ence between the diagonal components of the virial 
tensor Hat — IIt or, equivalently, an integral over 
the differential pressur e Pn — Pt 



7 = ^ (Hw - Ht) = J 



(PN - Pt) dy, (2) 



where 2A denotes the interfacial area of the two in- 
terfaces. The pressure calculation is based on the 
method proposed by Irving and Kirkwood,^^ but in 
contrast to their approach the pressure is not ho- 
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Table 1: Surface tension 7 in dependence of the slab thick- 
ness S (low temperatures). The total elongation of the sim- 
ulation box is indicated as £y in the direction perpendicular 
to the vapor-liquid interfaces and as £|| (= ix = £z) in the 
other spatial directions. The statistical error in terms of the 
final digit is shown in parentheses. 



mogeneously distributed between the positions of 
the two particles. To speed up the simulation the 
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Figure 4: Differential pressure pjy — Pt over the density for 
the temperature T = 0.7 and T = 1.0. The blue circles 
and brown triangles correspond to the minimum stable con- 
figuration for the corresponding temperature, while the red 
squares and green stars correspond to S = 12. The dotted 
line represents the zero line. 
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pressure is divided between the bins of the involved 
particles. It should be noted that this modification 
has a slight effect on the localized pressure tensor 
but leads to the same outcome for the overall sur- 
face tension.— 

Molecular dynamics (MD) simulations were con- 
ducted in the canonical ensemble with = 16 000 
particles. The equimolar thickness of the liquid slab 
5' was varied between 12 and the minimum stable 
configuration. The equimolar thickness was deter- 
mined using the saturated densities p' and p" for 
the given temperature, the simulation volume y, 
and the number of particles N 



i.e. S only depends on the boundary conditions ap- 
plied to the molecular simulation within the canoni- 
cal ensemble, not on the outcome of the simulation, 
and it does not vary over simulation time. The 
temperature was kept constant by a velocity scal- 
ing thermostat. All simulations were performed in 
a parallelepiped box. The elongation of the sim- 
ulation volume in y direction, i.e. normal to the 
interface, was ly = 50. 

For the simulation of a reference case the num- 
ber of particles N was increased to 300 000, the 
elongation in y direction was iy — 100 and a slab 
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Table 2: Surface tension 7 in dependence of the slab thick- 
ness S (high temperatures), cf. Tab. [Tl 
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Table 3: Differential pressure pjv — VT iu the center of the 
liquid slab in dependence on the slab thickness S (low tem- 
peratures), "n/a": data not determined; asterisks: liquid 
slab found to be unstable. 
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Table 4: Differential pressure pjv — Pt iu the center of the 
liquid slab (high temperatures), cf. Tab. [3] 



thickness of 5 = 40 was used. The temperature 
T was varied between 0.7 and 1.25, i.e. from the 
triple point temperature to 95% of the critical tem- 
perature. The simulations were carried out using 
the Isl molecular dynamics code.— The equation 
of motion was solved by a leapfrog integrator.— A 
time step of At = 0.002 was used. The equilibra- 
tion was conducted for at least 120 000 time steps, 
while the production ran for 840 000 time steps. 
The statistical errors given in the present study are 
equal to 3 times the standard deviation of 7 block 
averages, each over 120 000 time steps. 

A further series of simulations was conducted to 
validate whether relatively small cutoff radii are 
permissible when Janecek's cutoff correction is em- 
ployed. The simulation results support this conclu- 
sion, cf. Fig. [TJ 
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3. Results 

As described above, a series of simulations was car- 
ried out using a large liquid slab of S" = 40 in order 
to approximate bulk phase behavior. The resulting 
values 7oo and p'^ are used as a reference for the 
further simulations. The resulting surface tension 
is shown in Fig. [2] The regression 



7 = 2.94 



T 



1.23 



(4) 



is obtained, with Tc — 1.3126 according to Perez 
Pellitero et al.— The type of correlation is the same 
as proposed by Vrabec et al. for the truncated and 
shifted LJ potential, and their exponent (i.e. 1.21) 
is also very similar to the present one. 

Moreover, simulations were also performed for 
smaller slab thicknesses {S < 12). Thereby, MD 
runs were conducted with successively smaller val- 
ues of S, until a minimum stable thickness was 
reached for the given temperature. 

In Fig. [3J the density and the corresponding dif- 
ferential pressure profile for T = 0.7 is plotted over 
the y coordinate. It can be seen that the density 
in the center of the slab for 5 = 12 and 5 = 7 
almost matches the bulk liquid density at satura- 
tion, which is also plotted in Fig. [3] as reference. It 
is slightly smaller for S = 4.3. In the differential 
pressure the difference between the three simula- 
tions is more significant. The differential pressure 
in the center of the slab (y = 0) almost reaches the 
zero line for S = 12, while for 5 = 7 and S = 4.3 
the pressure tensor is anisotropic throughout the 
liquid slab. The differential pressure can be seen as 
an indicator for the fluid to be isotropic or homo- 
geneous, i.e. not influenced by the two interfaces. 
For S larger than 12, the differential pressure fluc- 
tuates around zero in the center of the liquid slab 
for T = 0.7. 

In Fig. 21 the plot of the differential pressure over 
the density is shown, which is obtained from the 
data displayed in Fig. [3] Additionally results for 
T = 1.0 are shown. The plot in Fig. |3] exhibits van 
der Waals loops in all cases. The red squares and 
green stars correspond to a large liquid slab, while 
the blue circles and brown triangles show the result 
of the smallest stable liquid slab. Like in Fig. |3l 
it is obvious that the differential pressure does not 
reach zero in the latter case. 

The resulting surface tensions are shown in Tabs. 
[T] and [5J For all temperatures the surface tension 
decreases when the liquid slab thickness decreases. 
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Figure 5: Reduced surface tension -ff-frxi over the slab thick- 
ness S for different temperatures. The dashed lines represent 
the expression 7/700 = 1 — / , where the temperature- 
dependent coefficients were adjusted to the simulation re- 
sults, yielding a(0.7) = 5.8, a(0.8) = 7.9, a(0.9) = 9.3, 
a(l.O) = 16, a(l.l) = 28, a{1.2) = 34, and a(1.25) = 93. 



The differential pressure in the center of the liq- 
uid slab, i.e. the minimum differential pressure, is 
shown in Tabs. [3] and SI As already discussed, the 
differential pressure in the center of the liquid slab 
increases with decreasing slab thickness. The den- 
sity corresponding to the minimum differential pres- 
sure, i.e. the maximum density, in the liquid slab is 
shown in Tabs. [S] and [HI At the largest slab thick- 
ness the value agrees with the bulk properties. For 
lower slab thicknesses the density does not reach the 
bulk values. The minimum stable slab thickness in- 
creases with increasing temperature. Reducing the 
slab thickness below this point results in a rupture 
of the liquid phase and a transition from planar to 
cylindrical or spherical symmetry. 

In Fig. \5\ the relative surface tension - reduced 
by 7oo(r) as obtained from the large slab simula- 
tions - is plotted over the slab thickness for differ- 
ent temperatures. Confinement between two planar 
vapor-liquid interfaces reduces the surface tension, 
and the numerical data suggest that this effect is of 
the order 1/5'^. In Fig. El the reduced density p' / 
p'oQ is plotted over the slab thickness for different 
temperatures. The relative density also decreases 
upon decreasing the slab thickness and, similar to 
the surface tension, this effect is approximately pro- 
portional to 1/5^ and becomes more significant at 
high temperatures. 

For a slab thickness S > 12, the surface tension 
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Figure 6: Reduced density p'/p'^ over the slab thickness S 
for different temperatures. The dashed hnes represent the 
expression p' / p'^ = 1 — b(T)/S^, where the temperature- 
dependent coefficients were adjusted to the simulation re- 
sults, yielding 6(0.7) = 3.0, 6(0.8) = 4.4, 6(0.9) = 6.3, 
6(f.O) = 12, 6(1. f) = 27, 6(1.2) = 37, and 6(1.25) = 85. 



agrees with the value for a large liquid slab, within 
the simulation uncertainty. For a fluid described by 
the LJ potential, e.g. methane,— this means that 
confinement effects are significant for slabs which 
are thinner than 4.5 nm. At high temperatures, 
the density in the center of the liquid slab deviates 
more significantly from the bulk value. 

The present results suggest that the reduction of 
the density and the surface tension due to confine- 
ment are related effects. In Fig. [3 the respective 
ratios are displayed together, which discloses an ap- 
proximately linear relationship. The regression 
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7oo 
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is obtained from the simulation results. 
4. Conclusion 

In the present work, molecular simulation was ap- 
plied to study the influence of the slab thickness 
on the interfacial properties for planar vapor-liquid 
interfaces. The present results prove that such an 
effect exists for thin slabs and quantifies it for the 
LJ fluid. The surface tension decreases with de- 
creasing slab thickness, and so does the density in 
the middle of the slab. The differential pressure 
does not reach zero for liquid slabs smaller than 12, 
which proves that under such conditions, a bulk- 
like region is absent. The confinement effects for 
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Figure 7: Reduced density p' /p'^ over the reduced surface 
tension 7/700 for different temperatures. The dash dotted 
line represents the regression from Eq. l(5j. 



the surface tension and the density were found to 
scale with in terms of the slab thickness S, so 
that a linear relation between both effects could be 
obtained. 

The present results depart from those obtained 
by Weng et al.— in a previous study, where no sys- 
tematic correlation between the slab thickness and 
the surface tension was found. For a LJ system at 
T = 0.818, Weng et al. detected minor fluctuations 
around a constant value (7 = 0.78 ± 0.02), without 
a clear tendency, for a range of slab thicknesses be- 
tween S = 5.0 and 9.0.— A juxtaposition with the 
present numerical data, cf. Tab. [T] and Fig. [5j ac- 
cording to which varying the slab thickness to such 
an extent has a significant influence on 7, clearly 
shows that there is a contradiction between present 
simulation results and the postulate of Weng et al. 
that "with film thickness . . . surface tension val- 
ues and density profiles show little variation."— For 
the simulations of Weng et al.,— however, no long- 
range cutoff correction was employed at all, and 
the computations were only carried out over 120 
000 time steps, as opposed to a million time steps 
for the present series of simulations. Since systems 
with an interface relax more slowly than the homo- 
geneous bulk fluid, the extremely short simulation 
time could consitute a serious limitation, affecting 
the accuracy of the results obtained by Weng et 
al.— to a significant extent. 

The comparison with results from a recent study 
of Malijevsky and Jackson— suggests that the 



6 



present results on confinenient by two parallel pla- 
nar vapor-liquid interfaces might also carry over 
qualitatively to confinement by the opposite sides of 
the single spherical interface that surrounds a small 
droplet. Therein, Malijevsky and Jackson come to 
the conclusion that for liquid drops, the size depen- 
dence of the surface tension is best described by 
two distinct, additive terms: The conventional Tol- 
man term, representing curvature, which increases 
the surface tension (i.e. the Tolman length is found 
to be negative), as well as "an additional curva- 
ture dependence of the l/R^ form" which causes an 
eventual decrease of the surface tension "for smaller 
drops."— Furthermore, Malijevsky and Jackson ob- 
serve that the characteristic droplet radius, below 
which this negative corrective term becomes domi- 
nant, "increases with increasing rc" and conjecture 
that "such a crossover occurs when ... no 'bulk' 
region can be assigned inside the drop. In this case 
even particles in the centre of the drop 'feel' the 
interface." 

The present results lend further plausibility to 
this conjecture of Malijevsky and Jackson. There 
could be a relation between their 1/R^ term and 
the confinement effect from the present study. 
According to such a hypothesis, these contributions 
would both represent the deviation from bulk-like 
behavior of the liquid phase due to confinement. 
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Table 5: Density p' in the center of the liquid slab in depen- 
dence on the slab thickness S (low temperatures). 
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0.49 (10) 


9 


0.62 (2) 
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Table 6: Density p' in the center of the liquid slab (high 
temperatures), cf. Tab.|5] 
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